Regular spiking in asymmetrically delay-coupled FitzHugh-Nagumo systems 
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Abstract — We study two delay-coupled FitzHugh- 
Nagumo systems, introducing a mismatch between the 
delay times, as the simplest representation of inter- 
acting neurons. We demonstrate that the presence of 
delays can cause periodic oscillations which coexist 
with a stable fixed point. Periodic solutions observed 
are of two types, which we refer to as a "long" and a 
"short" cycle, respectively. 

I. Introduction 

Being an inherent feature of a human, laziness was 
always that force which engendered invention of new 
devices, supposed to work instead of people. And 
in our epoch of vast technological progress, there are 
thousands of useful gadgets akeady existing. Though 
recently the science is advancing with seven-league 
strides, there are still a great number of phenomena 
which are waiting for a better insight. One has to 
admit that none of the existing complex machines 
and powerful computers can substitute a single human 
brain. Which means that we still do not draw close 
enough to clearing up a mystery of how this accumu- 
lation of grey matter really works. 

Since the end of the last century, study of neural 
networks picks up speed. In order to describe its in- 
tricate behavior, the brain is often represented as an 
ensemble of coupled nonlinear dynamical elements, 
capable of producing spikes and exchanging informa- 
tion between each other [1-3]. Such neural popula- 
tions are usually spatially localized and contain both 
excitatory and inhibitory neurons [4]. 

Some researchers, starting from the simplest case of 
two interconnected neurons, show how more compli- 
cated dynamics emerges in larger sets [5]. The others 
explore extremely complex network of subnetworks, 
focusing on the hierarchically clustered organization 



of interacting excitable elements [6]. 

Most studies base on the present oscillatory behav- 
ior of individual system elements, which then pro- 
duces observable patterns due to collective synchro- 
nization [7-9]. Thus, for modeling a single neuron, 
phase oscillators are often used. For instance, to char- 
acterize mutual dynamics of cells in certain brain ar- 
eas, responsible for giving the onset to Parkinson's 
disease or epilepsy, a well-known Kuramoto model is 
considered [10-13]. 

Here, we rely on the works by FitzHugh [14] and 
Nagumo et al. [15] who have shown that for describ- 
ing the main characteristics of a neuron dynamics, it 
is sufficient to consider a 2-dimensional system. The 
latter is also widely used nowadays as one of the sim- 
plest models for examining brain dynamics and has 
been essentially studied in many papers (see, for in- 
stance, [16-19] and references therein). 

Having an intention to move from simple to com- 
plex, we consider below a set of equations consisting 
only of two identical FitzHugh-Nagumo subsystems 
(see also [20,21]). Their interaction is described by 
a linear coupling term which includes delays (ri and 
T2), accounted for the fact that the signal transmission 
between neurons is not instantaneous: 
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Here {xi,yi) and {x2,y2) are the phase coordinates 
for the first and the second subsystem respectively. 
The parameter a determines whether the individual 
neuron is in the excitable regime or exhibits self- 



sustained periodic firing. Tiie time scale parameter e 
is ciiosen during the numerical simulations to be 0.01, 
which results in fast activator variables xi, X2 and 
slow inhibitor variables yi, y2. For further simplic- 
ity, the coupling strength C is also taken symmetric. 



II. Sketch dynamics 

A. Fixed point 

As it was already mentioned, the dynamics of an 
isolated 2-dimensional FitzHugh-Nagumo system is 
already well-investigated. Its single fixed point P2 = 
(-0,0^/3 —a) is stable for a > 1 and exhibits a super- 
critical Hopf bifurcation when the excitability param- 
eter crosses unity, which implies periodic spiking for 
o < 1. Provided that a > 1, the system is excitable, 
namely, if a sufficient external impulse is added, it 
emits a spike and then rests again in the P2 state. 

For our numerical simulations, we take a = 1.3, 
so that the individual subsystems are in the ex- 
citable regime. The coupling term of the consid- 
ered form is canceled for a fixed point orbit, thus, 
the 4-dimensional equilibrium P4 = (— a, — a + 
0^/3,— a,— o + a'^/3), being existent for the uncou- 
pled system, persists as well for the Eq. (dJ. Changing 
the coupling strength or the delays also does not influ- 
ence its stability, as it was recently shown [20]. 

B. Regular spiking 

However, besides the stable fixed point solution, 
the system (O can also produce periodic oscillations. 
Intuitively, this phenomenon can be explained as fol- 
lows. One can perturb, for instance, the first neuron, 
so that it emits a spike. Then, with the delay ri this 
perturbation reaches the second neuron, which pro- 
vokes it to spike as well. Again with the delay T2 the 
second neuron "informs" the first one that it has been 
stimulated, which causes a new run of the cycle, and 
the process repeats (see schematic representation in 
the Fig da)). 
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Fig. 1 

Schematic representation of periodic firing in 

A SYSTEM WITH DELAY. (A) "LONG" CYCLE; (B) 
"short" CYCLE 
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Fig. 2 

Phase portraits ((a) and (b)) and time series 
((c), (e) and (d), (f)) for long and short cycles 
respectively. On the time series plots, solid 
line correspond to xi^2 and dashed line 
correspond to yi^. the parameters are 
a = 1.3, £1,2 = 0.01, C = 0.5, n = 3, T2 = 1. 



Though, in the numerical simulations, starting from 
various initial conditions, we observed periodic solu- 
tions of two different types, which are referred to in 
the following as a "long" and a "short" cycle respec- 
tively. The former is of the period T ^ ti + T2, while 
the latter has the period T ^ (n + T2)/2. 

Again, intuitively, to obtain this second solution 
one would add an initial impulse not to one, but to 
both neurons, then roughly the short cycle dynamics 
can be plotted as in the Fig. [Hb). One could remark 
that, in this case, the initial perturbation for the second 
neuron should arrive before the delayed signal of the 
first one, namely for t € (0,ri). Although there are 
infinitely many variations for choosing the time mo- 
ment for the second impulse, in our numerical simula- 
tions we were able to observe only that pattern, which 
is depicted in the Fig.[Tlb). 

In the Fig. [2l we plot the phase portraits and the 
data series for these two attractors for certain fixed 
parameter values. 

III. Periodic solutions: deeper insight 

The next point to investigate in connection with the 
periodic firing patterns obtained, is a question whether 
these solutions exist for all couplings. Is their stability 



region large enough or such solutions appear only for 
separate parameter values? 

A. "Long " cycle 

As it was already noticed in [20], such oscillations 
appear through a saddle-node bifurcation of limit cy- 
cles, creating a pair of a stable and an unstable peri- 
odic orbit. In the Fig. Oa), the bifurcation curves of 
this attractor type are plotted in the (C, Ti)-plane, for 
T2 = 0.5, T2 = 1 and T2 = 2. It is easy to conclude, 
that with increasing T2 the bifurcation curve moves to 
the left, closer to the wall value C = 0. 




Fig. 3 

Bifurcation curves for appearance of the 
(n + t2)-peri0dic solution (long cycle) in 
(C, ri)-PLANE. (A) Diagrams for different t2 

VALUES, (b) Diagram for the case of n = t2 = t. 

(c) Overlay of the figures (a) and (b). Vertical 

DASHED line INDICATES A CRITICAL VALUE OF C. 



This implies that for some large enough coupling 
periodic firing still exists even if one of the delays is 
close to zero. For comparison, in the Fig.|3lb), the bi- 
furcation curve for the case ri = T2 is present. When 
overlaying the two graphs of (a) and (b) (Fig. He)), 
one can notice that the critical coupling value (indi- 
cated by a vertical dashed line) does not depend on 
the delay times difference, but only on their sum (see 
Appendix). 

In support to this last statement, we depict in the 
Fig. IHa, b) phase portraits and time series for three 
different periodic solutions, namely n = T2 = 2, 
Ti = 3,r2 = 1 and ti = 3.5, T2 = 0.5, while the 
sum of delays is always 4 and the coupling strength 
C = 0.5. As it could be clearly seen, the phase tra- 
jectories coincide perfectly as well as the time series. 
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Fig. 4 

Phase portrait and time series for 3 different 
long cycles with c = 0.5. 



We also would like to examine the question how 
the cycle period is related to the coupling terms. The 
Fig. |5la) represents several plots of the orbit period T 
vs. Ti, while T2 = 0.5, 1, 2 and C = 0.5. 

In the Fig. Ob), dependence of the period on C is 
depicted {t2 is the same as in (a), and ti is chosen so 
that the sum of delays does not change). As it is ex- 
pected (cf. [20]), T increases linearly with ri. How- 
ever, it decays with increasing C. 
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Fig. 5 

Evolution of the period for the long cycle 

SOLUTION DEPENDING ON THE COUPLING TERM. (A) 

Period T vs. ti for different fixed values of t2, 
C = 0.5. (b) Period T vs. C, for different values 

of n AND T2 SO THAT n + T2 = 4. 



B. "Short" cycle 

For the short cycle, the situation is almost the same. 
Again it is born through a saddle-node bifurcation. In 
the Fig. 13 a), we also plot the bifurcation curves, sepa- 
rating the regions of existence and absence of the short 
cycle, in the (C, ri)-plane (as earlier T2 = 0.5, r2 = 1 
and T2 = 2). Again with increasing T2 the bifurcation 
curve moves to the left, however, in comparison with 
the long cycle the short one occurs for larger values of 
coupling strength. And after laying over the curve for 
the case of equal delays ri = T2 (Fig. Ob)), one can 
notice that the critical C depends only on the delays 
sum (see Fig. Oc)). The phase portraits and the time 
series for three different periodic solutions, plotted in 
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Fig. 6 

Bifurcation curves for appearance of the 
(n + t2)/2-periodic solution (short cycle) in 
(C, ri)-PLANE. (A) Diagrams for different t2 
VALUES, (b) Diagram for the case of n = r2 = r. 
(c) Overlay of the figures (a) and (b). 



the Fig. Ul coincide as well. 
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Fig. 7 

Phase portrait and time series for 3 different 
short cycles with c 0.5. 



Finally, in the Fig. [8la),(b), the graphs disclosing 
the relation between the period and the coupling term 
configuration are presented. As in the case of the long 
cycle, T is a linear function of ri and has a gradual 
decrease on C. 

IV. Conclusions 

In the present paper we have considered two asym- 
metrically delay-coupled FitzHugh-Nagumo systems 
for modelling interacting excitable neural elements. 
Such an "intrusion" gives rise to the regular spiking 
in the system investigated. For sufficiently large cou- 
pling strength and delays, one can observe periodic 
solutions of two different types (long and short cy- 
cles), depending on whether only one subsystem is 
perturbed initially or both. The long cycle period ap- 
proximately equals n + T2, while the short one has a 



2 X, 4 0.5 c 1 

Fig. 8 

Evolution of the period for the short cycle 

SOLUTION depending ON THE COUPLING TERM. (A) 

Period T vs. n for different fixed values of t2, 
C = 0.5. (b) Period T vs. C, for different values 

OF Ti AND T2 so THAT n + T2 = 4. 



period of about a half of this amount. 

Furthermore, the numerical simulation, as well as 
the mathematical anlysis, shows that phase portraits 
and time series of these solutions do not depend on 
the difference of delays, but only on their sum. 
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Tliis corresponds to a system witli symmetric delay 
coupling, and the function X2{t) fully coincides with 
the function X2{t) of the initial problem, but with a 
shift along the time axis by At = (ri — T2)/2. 

We also note that the inhibitor variables yi, of 
the system ([T|l depend only on xi{t) and X2{t), respec- 
tively. Therefore, omitting them in the above analysis 
does not influence the resulting conclusion. 



Appendix 

a. transformation to symmetric coupling 
Consider the general system 

Xl= f{xi)+C{x2{t-T2)-Xl{t)), (2) 
X2 = f{002)+C{xi{t-Tl)-X2{t)). (3) 

Without losing generality assume that ti > T2 and 
denote t = {ti + T2)/2 and At = (ri — T2)/2, so 
that Ti = T + At and T2 = t — At. Then introducing 
a new function X2{t) = X2{t + At) we use 

X2{t - T2) = X2{t - t) 

in eq. ^ and rewrite the equation ^ as follows 

52(t) = f{x2{t)) + C{xi{t + At - Ti) - X2{t)) 
= f{x2{t))+C{xi{t-T)-X2{t)), 

which leads to 

xi = f{xi) + C{S:2{t-T)-xi{t)), 

(4) 

X2 = f{x2) + C{xi{t - t) - X2{t)). 



